
rm(list=ls())

library(tidyr)
library(tibble)
library(stargazer)
library(dplyr)
library(haven)
library(readstata13)
library(ggplot2)
library(cregg)
library(zoo)

#########################################################################
# REPLICATION FILES: MILITARIZATION AND PERCEPTIONS OF LAW ENFORCEMENT  #
# JOURNAL: BJPS                                                         #
# RESULTS IN SUPP. APP. FOR                                             #
# ROBUSTNESS CHECK FOR EFFECTIVENESS                                    #
# RESULTS AND FIGURE                                                    #          
#########################################################################


#### STEP 1: READ DATA                              ####

militconj <- read.dta13("MilitConjoint.dta")
baselines <- list()

#### STEP 2: RENAME VARIABLES                       ####

colnames(militconj)[colnames(militconj)=="uniform"] <- "uniform_numeric"
colnames(militconj)[colnames(militconj)=="weapon"] <- "weapon_numeric"
colnames(militconj)[colnames(militconj)=="gender"] <- "gender_numeric"
colnames(militconj)[colnames(militconj)=="race"] <- "race_numeric"

#### STEP 3: CHANGE TO FACTORS                      ####

militconj$Uniform<-factor(militconj$uniform_numeric,labels=c("Police", "Military"))
militconj$Weapon<-factor(militconj$weapon_numeric,labels=c("None", "Assault Rifle"))
militconj$Gender<-factor(militconj$gender_numeric,labels=c("Female", "Male"))
militconj$SkinColor<-factor(militconj$race_numeric,labels=c("Dark", "Light"))

#### STEP 4: RESULTS, FORCED CHOICE EFFECTIVENESS ####

cj(militconj, whicheffective ~ Weapon + Uniform +  Gender + SkinColor, id = ~idresp, estimate = "amce")

#### STEP 4: FIGURE, FORCED CHOICE EFFECTIVENESS ####

# UPLOAD RESULTS DATA
all_who <-  read.table("All_WhoEff.txt",sep="")

# DECLARE THEME
theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      axis.title.x =      element_text(size =12),
      legend.position = "none",
      panel.background = element_rect(fill = "#F5F5F5"),
      strip.background = element_rect(colour="black", fill="white")
    )
}

# FUNCTION TO PREPARE DATA FOR FIGURE
prepdata <- function(d){
  
  # prep estimates
  d$var <- rownames(d)
  colnames(d) <- c("pe","se","var")
  d$order <- 1:nrow(d)
  # compute Cis
  d$upper95 <-d$pe + 1.96*d$se
  d$lower95 <-d$pe - 1.96*d$se
  
  d$upper90 <-d$pe + 1.645*d$se
  d$lower90 <-d$pe - 1.645*d$se
  
  # define group
  d$group <- NA
  d$group[d$var %in% paste(c("1b",2),".FeatGender",sep="")]         <- "Gender"
  d$group[d$var %in% paste(c("1b",2),".FeatSkin",sep="")]           <- "Skin Color"
  d$group[d$var %in% paste(c("1b",2),".FeatUniform",sep="")]         <- "Uniform"
  d$group[d$var %in% paste(c("1b",2),".FeatWeapon",sep="")]       <- "Weapon"
  
  
  # order 
  d <- d[order(factor(d$group,levels=unique(d$group)[c(3,4,1,2)])),]
  d$order <- 1:nrow(d)
  
  # label attributes
  offset <- c("   ")
  
  d$var[d$group=="Gender"] <- paste(offset,c("Male","Female"))
  d$var[d$group=="Skin Color"] <- paste(offset,c("Light","Dark"))
  d$var[d$group=="Uniform"] <- paste(offset,c("Military","Police"))
  d$var[d$group=="Weapon"] <- paste(offset,c("Assault Rifle","None"))
  
  # sub in group labels
  dd <- data.frame(var= c("Uniform:",
                          " ",
                          "Weapon:",
                          "  ",
                          "Gender:",
                          "   ",
                          "Skin Color:"
  ),order=c(.5,2.1,2.5,4.1,4.5,6.1,6.5),
  pe=1,se=1,upper95=1,lower95=1, upper90=1, lower90=1,group=NA)
  d <- rbind(d,dd)
  d <-d[order(d$order),]
  d$var <- factor(d$var,levels=unique(d$var)[length(d$var):1])
  return(d)
}

# PREPARE DATA

all_who <- prepdata(all_who)
all_who$model<-"Who Effective"

# ORDER FACETS

ratings_all$model <- factor(ratings_all$model, levels = c("Effectiveness", "Civil Liberties",
                                                          "Corruption", "Neighborhood"))

# CREATE FIGURE A.2
yylab1  <- c("Effect on Pr(Security Personnel Selected as More Effective)")

all_whoeffective = ggplot(all_who, aes(y=pe,x=var)) +
  coord_flip(ylim = c(-.2, 0.2)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=4) +
  scale_y_continuous(name=yylab1) +
  scale_x_discrete(name="")

plot(all_whoeffective)

#pdf(paste("All_WhoEffective.pdf",sep=""),width=6,height=6) 
#all_whoeffective = all_whoeffective  + theme_bw1()
#print(all_whoeffective)
#dev.off()

#ggsave("All_WhoEffective.png", plot=all_whoeffective, width=6,height=6, units = "in", dpi = 300)

